Pressure correction in density-functional calculations 
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First-principles calculations based on density functional theory have been widely used in studies 
of the structural, thermoelastic, rheological, and electronic properties of earth-forming materials. 
The exchange-correlation term, however, is implemented based on various approximations, and this 
is believed to be the main reason for discrepancies between experiments and theoretical predic- 
tions. In this work, by using periclase MgO as a prototype system we examine the discrepancies in 
pressure and Kohn-Sham energy that are due to the choice of the exchange-correlation functional. 
For instance, we choose local density approximation and generalized gradient approximation. We 
perform extensive first-principles calculations at various temperatures and volumes and find that 
the exchange-correlation-based discrepancies in Kohn-Sham energy and pressure should be indepen- 
dent of temperature. This implies that the physical quantities, such as the equation of states, heat 
capacity, and the Griineisen parameter, estimated by a particular choice of exchange-correlation 
functional can easily be transformed into those estimated by another exchange-correlation func- 
tional. Our findings may be helpful in providing useful constraints on mineral properties at deep 
Earth thermodynamic conditions. 
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I. INTRODUCTION 

In recent decades, first-principles (FP) calculations 
with density functional theory (DFT) [U, Q have been 
widely adopted in studies of the structural, thermoe- 
lastic, rheological, and electronic properties of materi- 
als, molecules, and minerals in research areas that range 
from physics and chemistry to biology and geosciences 
S 11 M S 0, Hi- The modern implementation of FP 
techniques is based on different approximations, such as 
the pseudopotential (PP) method, the expansion of elec- 
tronic wavcfunctions by modal functions, and the choice 
of exchange-correlation (XC) functionals. As a result, 
the predictions made by different FP calculations may 
differ from one another moderately. 

Many techniques have been developed to enhance the 
accuracy and efficiency of FP calculations. Examples in- 
clude the projector-augmented wave (PAW) method JJ, 
and maximally localized Wannier functions [l0|. How- 
ever, due to inadequate knowledge, the exact form of 
exchange-correlation functional remains unknown and 
the search for a precise XC term has been an active area 
of research. The choice of XC in FP calculations has 
been shown to be the main reason for the inconsistencies 
between these various calculations pllIT^. [Tsl. [T^. [isl. [l^ . 
In particular, the pressures calculated by generalized gra- 
dient approximation (GGA) [l3] and local density ap- 
proximation (LDA) TS] are known to be systematically 
over- and underestimated, and a common way to choose 
the XC functional is to use the XC that gives the best 
agreement with experiments. 



In this article, we investigate the way in which such a 
pressure difference depends on volume and temperature. 
If this difference can be quantified, simulations can be 
performed using either exchange-correlation, and at the 
same time, obtaining the results of the other. This would 
be an efficient way of constraining the thermodynamic 
quantities. For example, Karki et al. pj| calculated the 
elastic properties of MgSi03 at pressures up to the lower 
mantle condition using LDA, which gave an equilibrium 
volume and elastic constant consistent with experiments 
[l9, 20]. However, calculations by Oganov et al. [T3| 
on MgSi03 have shown that the GGA gives an extremely 
accurate elastic constant and the volume dependencies of 
the elastic properties, although the pressure is overesti- 
mated. The authors thus applied a constant pressure 
shift to the equation of state (EOS) estimated by GGA 
and then matched the measured ambient pressure EOS, 
which resulted in an excellent match between the exper- 
imental observations and the FP calculation results. 



II. THEORY 

To quantify the relationship between the pressure dif- 
ference and the XC functionals, we first consider the total 
internal energy of a system, which is contributed by the 
interaction energy between particles and the total kinetic 
energy of the ions and electrons. At thermal equilibrium, 
this total internal energy is given by 
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where (• • • ) denotes the ensemble average, and {R/} and 
{ipn} represent ionic and electronic degrees of freedom. 
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Also, Ml and qi represent the mass and charge of the Ith 
ion. The last term i?Ks[{R-/}i {V'n}] is the Kohn-Sham 
energy functional of the ion-electron system. Within adi- 
abatic approximation, ii'KsilR-/}! {"0"}] is given by DFT, 
and i?KS depends on the configuration of ions for bulk 
systems. The first and the second terms in Eq. ^ are 
the classical kinetic and coulomb energy of the ions. For 
a crystal system with small oscillation of ions, one can 
write R/ = R5 + STLj, where and STLj are respec- 
tively the equilibrium position and the (small) displace- 
ment from Rj of the Ith ion. In addition, we can assume 
that each ion is under the influence of an effective local 
potential, given by {l/2)Ki^a{SRi,a)^- Here i^/,Q is the 
effective force constant, which includes the contributions 
from ion-ion and electron-ion interactions, and a is the 
component of displacement {a — 1, 2, 3). As a result, the 
total internal energy [Eq. ([T])] can be approximated by 
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which could be understood as the sum of zero- 
temperature energies (quantities with superscript 0) and 
the contribution due to thermal motions of ions (quanti- 
ties with (• • • )). 

In first-principles calculations, and {&E) depend 
on the choice of XC functional, whereas (T) equals 
3NkBT/2 regardless of the choice of XC, because the 
ionic motions are treated classically. Therefore, Eq. ([2]) 
can be written as. 



where 



U = -NknT+iE), 



{E) = Ut^ + E^s + m, 



(3) 



(4) 



is the total DFT energy of the system. We now consider 
the internal energy calculated by two different choices of 
XC, say XCi and XC2, the difference between U is given 
by 



AU{V,T) = Uxc^iV,T)~UxcAV,T) 



(5) 



As U^Q^ only depends on the equilibrium ionic configura- 
tion, it does not contribute to AU{V,T). Moreover, {6E) 
should be equal to 3NkBT/2, because SRj^a is regarded 
as a classical degree of freedom and 

(iif/,„((5i?,,„)2) = ifcsT, 

at thermal equilibrium. Therefore, AU becomes 

AUiV,T) = i?^s,xc, -£^kxc.: (6) 



which depends solely on the system volume. More impor- 
tantly, the difference between the calculated pressures, 
given by 
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(where S represents the entropy of the system) should be, 
in general, a function of V and T. However, according to 
Eqs. ([5]) and although both -ExCi ^-nd Exc2 depend 
on temperature, their difference, Exci ~Exc2 should not. 
Consequently, the resultant difference in pressure should 
also be temperature independent, that is. 



Pxc^{V,T) - Pxc2{V,T) = APi2{V). 



III. METHODOLOGY 



(8) 



To illustrate the relationships in Eqs. © and ([5]), we 
perform FP calculations on periclase MgO at thermo- 
dynamic conditions that range from ambient conditions 
to pressures and temperatures that are close to those of 
a deep Earth environment. Periclase MgO has a face- 
centered cubic structure and is known to be stable in 
the pressures and temperatures being considered in this 
work. Therefore, structural changes, such as phase tran- 
sition and melting, are avoided. 

The EOS at each temperature are determined by first- 
principles molecular dynamics simulations. For each sim- 
ulation, the simulation supercell consists of 64 atoms (32 
MgO units). These atoms were first arranged in a face- 
centered cubic (fee) configuration, and the atomic tra- 
jectories and electronic orbitals were evolved via Car- 
Parrinello molecular dynamics (CPMD) [2l|, with the 
temperature controlled by the Nose thermostat technique 
at T = 300 K, 1000 K, 2000 K,... 4000 K, at var- 
ious cell volumes. The calculations are repeated using 
both LDA and GGA exchange-correlations. The pres- 
sures determined at each temperature are fitted against 
the third-order Birch-Murnaghan EOS ,23] ■ As a result, 
two EOSs at each temperature are obtained, one LDA 
and one GGA. 

The choice of the XC energy functional is implemented 
in the generation of the pseudopotentials. In this study, 
two pairs of pseudopotentials for Mg and O are gener- 
ated, one pair generated using LDA and the other us- 
ing GGA. Apart from this difference in XC, all of the 
other parameters, such as the cutoff radius and the va- 
lence states, are the same. For Mg, a norm conser ving 
pseudopotential 24] with nonlinear core correction [29] 
is used, whereas for O an ultrasoft [2^ pseudopotential 
is used. The reference configuration is 3sl.5 3p0 3d0 
for Mg and 2s2 2p2 for O. For CPMD, a fictitious elec- 
tron mass /ig — 400 rUe and a simulation time step of 
At = 12 atomic time units 0.3 fs) are used. Elec- 
tronic wavefunctions are expanded by planewaves with a 
cutoff of 30 Ry, and the cutoff for charge density is 240 
Ry. To ensure enough statistical data, the simulations 
are run for more than 4 ps at each V , T point. 
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IV. RESULTS AND DISCUSSION 
A. Equation of states and bulk modulus 

The EOS given by our LDA and GGA calculations 
are shown in Fig. 1. The results are fitted against the 
third-order Birch-Murnaghan EOS, with Vq = 3782.30 
a|, = 177.486 GPa, and K'^ = 4.026 for LDA and 
= 149.320 GPa, and K'^ = 4.080 



= 4046.43 a%, 



for GGA. To ensure that our calculations are compatible 
with existing results, we compare our results with those 
obtained by Karki et al. [131, and Isaak et al. [11]. The 
bulk moduli for each EOS are shown in Fig. 2. We re- 
calculate the static EOS of Karki et al. [l^l by using the 
same PP files as in their work [2^, and the recalculated 
static EOS is in excellent agreement with their published 
EOS parameters. In addition, such EOS was shown to 
be highly consistent with the experimental measurements 
at moderate pressures (P < 170 GPa). Therefore, Fig. 1 
and Fig. 2 illustrate the discrepancy that arises from dif- 
ferent methodologies, such as the implementation of PP, 
the choice of XC functionals, and the interaction poten- 
tials. 
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FIG. 1: (Color online) Calculated equations of state for MgO 
by LDA (circle) and GGA (star) at K. The results of Karki 
et al. [23] (dashed line) and Isaak et al. [2^ (dotted line) 
are given for comparison. The EOSs are drawn using the 
published EOS parameters. 

At first glance, this discrepancy appears to reduce as 
the system volume increases. This is because both the 
pseudopotential errors and the difference between the XC 
functionals become less significant as atomic separation 
increases. Also, the GGA-calculated pressure and the 
bulk modulus in this work deviate systematically from 
those of LDA. In addition, our GGA results are more 
compatible with those estimated in Ref. 27 and Ref. 28 
than our LDA calculations. This indicates that the GGA 
calculation should, in principle, be closer to the experi- 
mental measurements. The Mg PP used in Ref. 27 was 
optimized by using a combination of various electronic 
configurations, and a norm-conserving O PP was used. 
Nevertheless, the optimization of a PP configuration is 
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FIG. 2: (Color online) Isothermal bulk moduli given by vari- 
ous calculations. 
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FIG. 3: (Color online) DFT total energy difference between 
GGA and LDA against volume at different temperatures. The 
results for T = K are joined for visualization. 



the subject of further study. We concentrate on the er- 
rors introduced by the XC alone. 



B. Energy and pressure differences 

According to Eq. ([5]), it is intuitive to calculate the 
difference in energy at each temperature and volume. In 
Fig. 3, we plot the difference in energy, IS.E = Eqqa — 
-EldAj against V at each temperature and find that it 
does not depend on temperature. This is consistent with 
the prediction of Eq. The differences between the 
pressures estimated by LDA and GGA at each temper- 
ature, AP — Pgga — Plda, are shown in Fig. 4. As 
can be seen, the differences in pressure at all temper- 
atures almost coincide, with a maximum deviation of 
about 1 GPa, which is much smaller than the typical 
statistical error in pressure estimated by the DFT-based 
calculations. This is again in accordance with the con- 
jecture of Eq. ([H]), and the pressure difference should be 
independent of temperature. Within the pressure range 
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FIG. 4; (Color online) Difference in pressure AP = Pqga ~ 
Plda at each temperature. Note that the pressure differ- 
ence is in the order of 10 GPa, and the maximum difference 
between these points is about 1 GPa, which is negligible in 
DFT-based calculations. The pressure difference for T — K 
is fitted against AP = ai/V + a2/V'^ with ai = 1.53 x 10* 
GPa a% and 02 — 1.06 x 10* GPa a% for visualization. 



we have examined (0 GPa < P < 170 GPa), AP{V) 
behaves asymptotically as a polynomial of 1/V, that is, 



AP{V) 



(9) 



Here, a few comments are in order. The present anal- 
ysis implies that the thermodynamic quantities, such as 
internal energy and pressure, calculated by FP method- 
ologies with a particular XG functional can easily be 
transformed into those obtained by another choice of XG 
functional. For example, the difference in heat capacity, 
given by 



ACv = Cy^GGA — Cy.LDA 

d{EGGA - ^lda) d{AE) 



dT 



dT 



0, (10) 



at each volume should vanish because AE does not de- 
pend on temperature. This implies that the heat capacity 
calculated by any XG should be the same. In addition, 
the Griineisen parameter. 
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calculated by different XG functionals should also be 
the same. Here a = {l/V){dV / dT)p and Kt = 
—V{dP/dV)T are, respectively, the coefficient of thermal 
expansion and the isothermal bulk modulus. To prove 
the statement, we consider the following thermodynamic 
relation. 



■^^ .^m (s ^«A-. ,12, 



dT 



dT J p\dVy J, 



Therefore, Eq. (Ilip can be written as 

As a result, the difference in the Griineisen parameter. 
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should vanish either. 

The present estimation of pressure and energy differ- 
ences is based on calculations of the cubic MgO structure. 
In principle, our analysis should be valid for any crystal 
system with small oscillation of ions, although the corre- 
sponding FP calculations may require extra care in the 
estimation of pressure. Although not shown in this work, 
our preliminary calculations on MgSiOa perovskite and 
post-perovskite phases [sOl support our proposed theory. 
For a non-cubic crystal structure, a finite strain may in- 
troduce shear stress to the system. In such a case, the 
EOS should be determined by cell dynamics algorithms 
[3, It should also be noted that our conclusion on 
the pressure correction relies on the assumption that the 
thermal contribution to the total internal energy is gov- 
erned by classical mechanics. However, this may not be 
valid in situations where quantum mechanical interac- 
tions, such as phonon vibration should be taken into ac- 
count. In particular, at extremely low temperatures, long 
range acoustic phonon modes play an important role in 
various thermodynamic phenomena. Nevertheless, recent 
works on structural phase transition in MgSiOa at core- 
mantle boundary conditions [sij and postspinel transi- 
tion in Mg2Si04 32] have shown that high-temperature 
thermodynamic properties estimated by phonon-based 
calculations should not depend on the choice of XG func- 
tionals, as long as quasiharmonic approximation (QHA) 
remains valid in the thermodynamic conditions of inter- 
est. For example, different choices of XG functionals only 
affect the position of the estimated phase boundaries in 
these works but not the Glapeyron slopes. In addition, 
the recent theoretical study of ultrahigh-pressure EOS of 
MgO [3^ has also shown explicitly that, within the QHA 
validity regime, the difference between a LDA-EOS and a 
GGA-EOS is independent of temperature. More impor- 
tantly, it should be noted that QHA requires the vibra- 
tional amplitude of each atom to be small, which is the 
major assumption of our pressure correction conjecture. 
As a result, the work of Wu et al. [sll ubiquitously sup- 
ports our hypothesis of pressure correction, even the FP 
methodologies used (lattice dynamics) in this work are 
different from ours (Gar-Parrinello molecular dynamics). 

The validity and behavior of AP{V) are also subjects 
for further investigation. In the pressure range we have 



5 



investigated, AP{V) is finite and is well approximated 
by Eq. which implies that the pressure difference 
should increase with a decreasing volume. However, this 
is not valid either when V approaches infinity (V oo) 
or when V becomes vanishingly small {V ~ 0). In the 
former case, the system consists of isolated atoms, and 
the energy difference that is due to different XC is con- 
stant; as a result AP{V) is zero as y ^ oo. In the latter 
case, the inter-atomic distance approaches zero, the in- 
teracting electrons should be well described by the free 
electron gas model, and LDA and GGA should give the 
same pressure. Last but not the least, the present study 
focuses on the pressure correction due to the chosen XC 
functional, which is in opposite to the idea of volume cor- 
rection suggested by Wu et al. [ss'l, in which an EOS is 
corrected to match experimental data. A combination of 
these EOS corrections will make various FP-based results 
more transferable, thus allowing one to have an effective 
recipe to transform the calculated thermodynamic quan- 
tities from one condition to another. 



V. CONCLUSION 

To conclude, by using MgO as a prototype system, we 
have studied the discrepancy in pressure that is due to 



different choices of XC functionals in DPT calculations. 
We have found that the differences in energy and pressure 
for GGA and LDA calculations should be independent of 
temperature. As a result, one may easily constrain the 
XC error at arbitrary temperatures. This may lead to 
a better estimation of thermodynamic properties, such 
as heat capacity and the Griineisen parameter, in the 
systems that are being investigated by the PP and high- 
pressure communities. 
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